!> author: 左志华
!> date: 2022/05/14
!>
!> 卡方检验
!>
!> 网址：https://rosettacode.org/wiki/Verify_distribution_uniformity/Chi-squared_test#Fortran

module gsl_mini_bind_m
 
    use iso_c_binding
    implicit none
    private
 
    public :: p_value
 
    interface
        function gsl_cdf_chisq_q(x, nu) bind(c, name='gsl_cdf_chisq_Q')
            import
            real(c_double), value :: x
            real(c_double), value :: nu
            real(c_double) :: gsl_cdf_chisq_q
        end function gsl_cdf_chisq_q
    end interface
 
contains
 
    !> Get p-value from chi-square distribution
    real function p_value(x, df)
        real, intent(in) :: x
        integer, intent(in) :: df
 
        p_value = real(gsl_cdf_chisq_q(real(x, c_double), real(df, c_double)))
 
    end function p_value
 
end module gsl_mini_bind_m

program chi2test
 
    use gsl_mini_bind_m, only: p_value
    implicit none
 
    real :: dset1(5) = [199809., 200665., 199607., 200270., 199649.]
    real :: dset2(5) = [522573., 244456., 139979., 71531., 21461.]
 
    real :: dist, prob
    integer :: dof
 
    write (*, '(A)', advance='no') "Dataset 1:"
    write (*, '(5(F12.4,:,1x))') dset1
 
    dist = chisq(dset1)
    dof = size(dset1) - 1
    write (*, '(A,I4,A,F12.4)') 'dof: ', dof, '   chisq: ', dist
    prob = p_value(dist, dof)
    write (*, '(A,F12.4)') 'probability: ', prob
    write (*, '(A,L)') 'uniform? ', prob > 0.05
 
    ! Lazy copy/past :|
    write (*, '(/A)', advance='no') "Dataset 2:"
    write (*, '(5(F12.4,:,1x))') dset2
 
    dist = chisq(dset2)
    dof = size(dset2) - 1
    write (*, '(A,I4,A,F12.4)') 'dof: ', dof, '   chisq: ', dist
    prob = p_value(dist, dof)
    write (*, '(A,F12.4)') 'probability: ', prob
    write (*, '(A,L)') 'uniform? ', prob > 0.05
 
contains
 
    !> Get chi-square value for a set of data `ds`
    real function chisq(ds)
        real, intent(in) :: ds(:)
 
        real :: expected, summa
 
        expected = sum(ds)/size(ds)
        summa = sum((ds - expected)**2)
        chisq = summa/expected
 
    end function chisq
 
end program chi2test

! Dataset 1: 199809.0000  200665.0000  199607.0000  200270.0000  199649.0000
! dof:    4   chisq:       4.1463
! probability:       0.3866
! uniform? T
!  
! Dataset 2: 522573.0000  244456.0000  139979.0000   71531.0000   21461.0000
! dof:    4   chisq:  790063.2500
! probability:       0.0000
! uniform? F